#load MLM libraries
library(nlme)
library(lme4)

#Global settings
setwd("./")

#Reading in the Data
legislator.data <- read.delim("TobaccoLegislators.dat",sep="\t",dec=".",header=T)
legislator.data$legislator <- 1:527
attach(legislator.data)

##Center predictors
house <- ifelse(house==0,-1,house) #Senate was coded with a "0" in raw dataset
money <- money - mean(money)
acres <- acres - mean(acres)

##Aggregation Demo
#Regression: Moderated Regression
moderated.regression <- lm(voting~money*acres)
plot(effect("money*acres",moderated.regression))
summary(moderated.regression)

#Accounting for Hierarchy: Aggregation.
AGGREGATE <- data.frame(state = (levels(state)))
AGGREGATE$money <- tapply(money,state,mean)
AGGREGATE$acres <- tapply(acres,state,mean)
AGGREGATE$voting <- tapply(voting,state,mean)
aggregate.regression <- lm(AGGREGATE$voting~AGGREGATE$money*AGGREGATE$acres)
summary(aggregate.regression)

#MULTILEVEL MODELLING!
#Multilevel model with random intercept
mixed <- lme(voting~money,random=~1|state); summary(mixed) #Using the lme() function
mixed <- lmer(voting~(1|state)+money); summary(mixed) #Using the lmer() function

#Moderated multilevel model with random intercept and random slope
mixed <- lme(voting~money*acres,random=~1+money|state); summary(mixed) #Using the lme() function
mixed <- lmer(voting~(1|state)+money*acres+(money|state)); summary(mixed) #Using the lmer() function

#Moderated multilevel model with random intercept and a covariate
mixed <- lme(voting~house+money*acres,random=~1|state); summary(mixed)  #Using the lme() function
mixed <- lmer(voting~(1|state)+house+money*acres); summary(mixed)  #Using the lmer() function

#Simple effects testing in multilevel models
low.acres <- acres + sd(acres)
high.acres <- acres - sd(acres)
low.money <- money + sd(money)
high.money <- money - sd(money)

#Simple effects testing
mixed.low.acres <- lme(voting~house+money*low.acres,random=~1|state); summary(mixed.low.acres)
mixed.high.acres <- lme(voting~house+money*high.acres,random=~1|state); summary(mixed.high.acres)
mixed.low.money <- lme(voting~house+low.money*acres,random=~1|state); summary(mixed.low.money)
mixed.high.money <- lme(voting~house+high.money*acres,random=~1|state); summary(mixed.high.money)

#Baseline Model to calculate Effect Size.
mixed.baseline <- lme(voting~1,random=~1|state); summary(mixed.baseline) #Using the lme() function
mixed.baseline <- lmer(voting~(1|state)); summary(mixed.baseline) #Using the lmer() function

##### Advanced Applications #####
#3-level model with states nested within parties
party <- as.factor(party)
mixed <- lme(voting~house+money*acres,random=~1|party/state); summary(mixed) #Using the lme() function
mixed <- lmer(voting~(1|party)+(1|state:party)+house+money*acres); summary(mixed) #Using the lmer() function

#Nested Growth Curve Models - (HYPOTHETICAL DATA) Pretend that each legislator's vote on all pro-tobacco legislation was recorded (instead of the overall percentage) and the sequential order of each vote was stored in a variable called "time".
mixed.growth <- lme(voting~money*time,random=~1+time|state/legislator) #Using the lme() function
mixed.growth <- lmer(voting~(1|state)+(1|state:legislator)+money*time*(time|state)*(time|state:legislator)) #Using the lmer() function

#Generalized Multilevel Models - Poisson
mixed.poisson <- lmer(n_reps~(1|state)+gsp,family=poisson(link=log));summary(mixed.poisson)

#Generalized Linear Mixed Models - Multilevel logistic regression
mixed.bernoulli <- lmer(party~(1|state)+voting,family=binomial(link=logit));summary(mixed.bernoulli)

#Multilevel model with cross-classification
mixed.cross.class2 <- lmer(voting~(1|party)+(1|state)+house+money*acres); summary(mixed.cross.class2) #Cross-classification in 2 groups
mixed.cross.class3 <- lmer(voting~(1|party)+(1|state)+(1|house)+money*acres); summary(mixed.cross.class3) #Cross-classification in 3 groups
